Adaptive Lanczos-vector method for dynamic properties within the DMRG 
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Current widely-used approaches to calculate spectral functions using the density-matrix renormal- 
ization group in frequency space either necessarily include an artificial broadening (correction-vector 
method), have limited resolution (time-domain density-matrix renormalization group with Fourier 
transform method), or are limited to low-energy properties or single dominant modes (original con- 
tinued fraction method). Here we propose an adaptive Lanczos-vector method to calculate the 
coefficients of a continued fraction expansion of the spectral function iteratively. We show that one 
can obtain a very accurate representation of the spectral function very efficiently, and that one can 
also directly extract the spectral weights and poles for the discrete system. 
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Dynamical quantities such as the local density of 
states, the single-particle spectral weight, or the dynami- 
cal spin or charge correlation functions are of central im- 
portance in theoretical and experimental condensed mat- 
ter physics. Since electrons in solids are interacting quan- 
tum objects, a reliable calculation of their properties usu- 
ally has to resort to numerical approaches. The density- 
matrix renormalization group (DMRG)i is one such algo- 
rithm. Within the DMRG, the calculation of dynamical 
quantities is a considerable challenge. A first attempt 
by Hallberg was based on a continued fraction expansion 
(CFE)i^ Subsequently, Kiihner and White showed that 
this approach is suitable, "if only the low-energy part of 
the correlation function is of interest, or if the bulk of 
the weight is in one single peak" ^ and applied the cor- 
rection vector method^i^ which since then has success- 
fully been applied to many model systemsi^i^ However, 
this approach has the drawback that one needs to in- 
troduce an artificial broadening into the spectra, which 
can be viewed as convolution of the true spectral func- 
tion with a Lorentzian of width 77. As one is eventu- 
ally interested in the limit 77 — >■ 0, one has to "deconvo- 
lute" the spectrum at the end of the calculation. This is, 
like analytic continuation of Monte-Carlo data, a numeri- 
cally ill-defined procedure.— Furthermore, the calculation 
of the correction vector at every step of the DMRG is 
very time-consuming and the frequencies that one can 
address are limited. Another popular approach to cal- 
culating spectral functions is to Fourier-transform time- 
dependent DMRG datafi To obtain good frequency res- 
olution, one has to calculate time-domain data over a 
long time interval. However, accessing long time scales is 
limited by either a loss of accuracy due to the approxi- 
mate nature of the DMRG, or by finite-size effects such as 
reflections from open ends. A method that can resolve 
spectral features with high resolution and no artificial 
broadening is therefore desirable. Here we present an 
adaptive Lanczos-vector method (ALM) that takes ad- 
vantage of the CFE used in Rcf. [2, but gains efficiency 
by adapting the basis as in adaptive time evolution. 



The spectral function for the operators A and B is 
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G A b(^) denotes the zero-temperature Green's function. 
Here we will take B ^ A' , yielding 
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where lipo) is the ground state of Hamiltonian H, Eq is 
the ground-state energy, and s = +1/ — 1 when A is 
a bosonic/fermionic operator. One of several ways to 
represent the resolvent in Eq. ^ is the CFE^i^ 
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with a similar expression for G^^^ . The coefficients , bi 
of this CFE can be calculated using the recursion formula 
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which is essentially the one used in the Lanczos method, 
but with starting vector |/o) replacing a random vector. 
Thus, we will call the Lanczos vectors in the fol- 
lowing. Hallberg^ used the CFE ^ and the recursion 
formulae ^ to obtain spectral functions via a multi- 
target DMRG, i.e., by optimizing the DMRG basis for 
the ground state and the Lanczos states simultaneously. 

The structure of the recursion formula ([3]) suggests im- 
plementing an iterative method to calculate the Lanczos 
vectors in the DMRG that optimizes the basis for only the 
three Lanczos vectors needed at each recursion step. One 
initially calculates the ground state ji/'o) of the system to 
the desired accuracy with the usual DMRG algorithm. 
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One then performs single finite-system sweeps, simulta- 
neously targeting the additional vectors |/o) and 
from which one can evaluate oq, ai, and &i. The recursion 
proceeds by replacing \tpo), |/o), and |/i) by |/o), 
and the new vector I/2}. At this point, a technical sub- 
tlety arises: |/o) and cannot be recalculated because 
there is no longer a condition to optimize them. Instead, 
we transform the wave function from the previous finite- 
size DMRG step to the new superblock configuration^ 
at every step of the DMRG sweep. One DMRG sweep 
through the system thus suffices to calculate the Lanczos 
vector l/j) and the parameters ai and bi. As one needs 
to target only three Lanczos vectors simultaneously, the 
number m of basis states necessary to obtain an accurate 
representation is generally substantially smaller than in 
the original algorithmi^ We will come back to this point 
later. Here we emphasize only that we avoid calculating 
the ground state at each DMRG step, speeding up the 
calculation dramatically; the most time-consuming part 
left in the Lanczos iteration is now the diagonalization of 
the reduced density matrix needed to initialize the next 
DMRG step. After iterating the recursion relation often 
enough, one obtains a sequence {(fli, from which one 
can calculate the Green's function, Eq. ([2]). 

We test the method on a model of spinless fermions on 
a chain with Hamiltonian 

H = -J2 ('=1+1 + + ' 

i i 

where c^-^^ denotes the usual creation (annihilation) op- 
erators for an electron at site i, and rii denotes the oc- 
cupation number operator. Here we take A = Ci and 
calculate the spectral function p{x = i,uj) at the center 
of the chain for a system at half filling. In Fig. [TJ we 
compare the original (OLM) and the ALM implementa- 
tion of the CFE with the exact solution on a chain of 
length i = 40 in the free-fermion limit, [/ = 0. Both 
implementations do not give a clear hint as to when to 
stop the Lanczos iteration. We have evaluated the broad- 
ened spectral functions after 200 Lanczos iterations, see 
below. Following Kiihner and White,— we do not target 
all Lanczos vectors in the OLM but only a small fraction 
of them. We then calculate the remaining vectors by 
straightforward application of the recursion. To achieve 
good convergence of the DMRG, we assign a weight of 0.5 
to the ground state and equal weights to the other Lanc- 
zos vectors (one could also give them weights according 
to their spectral weights). This large freedom, especially 
in choosing the number of targeted Lanczos states, is in 
our a view a disadvantage of the OLM, as the optimal 
parameters are not obvious, and the results depend on 
their choice. Employing the truncation error as measure 
of the quality of the spectral function is also not possible, 
as it depends only on the number of targeted states and 
their weights. 

In Fig.[T]we have used m — 1500, 25 sweeps, and eleven 
target states, yielding a maximum truncation error of 
10~^. There is no way to determine the quality of the 
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FIG. 1. (Color online) Broadened spectral function for spin- 
less fermions (L = 40, f/ = 0, 77 = 0.05, 200 Lanczos vectors) 
calculated with the new adaptive implementation (ALM), the 
original implementation (OLM) and compared to the exact 
solution. The DMRG truncation number is given by m. 



representation of the following Lanczos states (one can 
only check their orthogonality), but if we would have 
tried to target all Lanczos states, the truncated weight 
would have been be much higher than 10~^. In addition, 
the convergence is unstable, as small changes in the first 
Lanczos vectors will lead to bigger changes in subsequent 
Lanczos vectors. 

In the ALM, we can increase the number of states dur- 
ing the Lanczos iteration, as we do not need to calculate 
the ground state at each step. We notice that the num- 
ber of states needed for a good representation of Lanczos 
vectors increases with the number of iterations desired. 
For m = 1500, the maximum truncation error per sweep 
for the ALM increases from 10~^^ (ground state and first 
two Lanczos vectors) to 10^** for the last three Lanczos 
vectors. As already noted for the OLMjS, the low-energy 
portion (w ^ 1) is well-reproduced by both methods. De- 
viations occur in the high-energy part of the spectrum, 
notably in both position and weight distribution for the 
OLM. Here the ALM already shows much better agree- 
ment with the exact solution for m = 1500. Increasing to 
m = 3000 for the Lanczos iterations in the ALM, one al- 
most perfectly reproduces the exact solution. Note that 
m = 1500 is the maximum number of states accessible in 
the OLM with our computational resources; the simul- 
taneous optimization of the ground state and all other 
Lanczos states prevents calculations with larger m. 

Another advantage of the ALM is its shorter run-time. 
For TO = 1500 our calculations used ~ 30ft. on a standard 
workstation for the ALM and « 90/i for the OLM. We 
emphasize, however, that the run-time is determined by 
several parameters and thus should be interpreted with 
some care. For example, the ALM scales linearly with 
the number of Lanczos iterations, whereas the original 
implementation is roughly independent of the iteration 
number. The run-time of the latter, however, strongly 
depends on the number of Lanczos states targeted, be- 
cause more target states increase the number of DMRG 
sweeps needed to achieve convergence. 

Since the CFE gives an analytical expression for the 
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FIG. 2. (Color online) Spectral weights evaluated for difTerent 
numbers of coefficients (a) and flow of the eigenvalues as a 
function of the number of Lanczos steps i (b) for spinless 
fermions {L = 40, U = 0, m = 3000). 



spectral function, wc can take the limit 77 — > so that 
the spectral function p^^^{bj) = J2n^n ^('^ ~ ^n) — 
s J2n <5(w + tj-'n) Is & serics of (5- functions, where a;„ = 
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states of H with cigcnenergies i?„ . Evidently, all of the 
information on the spectral function is contained in the 
weights 51^ and poles u>n = ±Ai?„. Furthermore, one 
will always obtain a discrete and finite set of poles for 
any finite system. Since the DMRG treats finite systems, 
the spectral weights and pole positions can be calculated 
directly. Such a calculation has two important advan- 
tages: (i) One can study the size-dependence of spectral 
properties in a very controlled manner, (ii) It is possible 
to obtain precise values for the energies of the low-lying 
excitationsjiii they are given by the pole positions a;„, 
which can be directly obtained as the eigenvalues of the 
tridiagonal matrix Tij = ai-iSij + hi ((5j,i+i -f (5i+ij)4 
The poles of the Green's function are located on the real 
axis and they are discrete. Therefore, one can integrate 
along a closed path C in the complex plane chosen to en- 
close one single singularity, say at uin- According to the 
residue theorem, the weight f2„ of this pole is then 
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[G{uin - e + 17) - G{ujn + e + 17)] d'j . 



The parameter e must be chosen to be smaller than the 
distance to the next eigenvalue, e < |w„ — a;„±i|. Using 



this procedure with the position of the poles known, the 
spectral weights can be calculated to high precision using 
numerical integration. 

In Fig. HJa) we compare the weights and positions of 
the poles calculated with the DMRG and the ALM to 
the exact values on a chain of length i = 40, which are 
given by fi„ = ;^|sin((fc/ -I- fc„)a;)|^, iOn = 2cos(fc„), 
where fc„ = -^fi, n e {!,..., L}, and kj is the Fermi 
wave vector. For small energies, the agreement is nearly 
perfect, while deviations occur for large energies. The 
quality of the agreement depends on the number of co- 
efficients a.i, bi taken into account. One finds that after 
an initial improvement, no improvement occurs when one 
further increases the number of coefficients. 

The origin of this behavior can be understood from 
Fig. [DJb), which depicts the fiow of the eigenvalues of 
the matrix Tij as a function of the number of Lanczos 
coefficients. The convergence of the first few eigenvalues 
is evidently rapid. The flow of the eigenvalues therefore 
gives a nice criterion for stopping the iterations (also for 
the OLM). However, after approximately 50 iterations, 
an eigenvalue with nearly vanishing weight appears and 
subsequently moves rapidly to zero energy. More such 
eigenvalues follow, at rapidly increasing frequency. This 
appearance of so-called "ghost" eigenvalues is well-known 
in the Lanczos method. These ghost eigenvalues are 
caused by the loss of orthogonality of the Lanczos vec- 
tors due to numerical errorJii^, Their appearance here 
is therefore not surprising, in particular because the cal- 
culated Lanczos vectors also include an error from the 
approximate DMRG representation. We find that this 
effect is enhanced as one reduces the number of states 
within the DMRG. Thus, the study of the fiow of the 
eigenvalues makes it possible to control the quality of 
the spectral function. Note that such ghost eigenvalues 
also occur within the OLM. As far as we know, this prob- 
lem was never addressed in detail for this method. While 
the ghost eigenvalues seem, at first glance, to be a seri- 
ous problem, we emphasize that, for the spectral func- 
tion in particular, they do not appear to cause real harm 
because they possess only very small weight. This is evi- 
dent from the scaling in Fig. [Hb) and also from the fact 
that all the ghost eigenvalues are located on the abscissa 
in Fig. [2l^a). This observation is just an empirical one 
at present. However, as long as this remains true, the 
method will be insensitive with respect to the occurrence 
of ghost eigenvalues, except for regions with very small 
spectral weight. Here it is difficult to distinguish ghosts 
from small but real spectral weights. One possibility is 
to examine their convergence Jii^ Another problem with 
ghosts is that they lead to a violation of sum rules due 
to double-counting. If one adds up just the real spectral 
values there will be a missing weight that is given by the 
sum of all ghost values. As long as the weight of the ghost 
values is small, this will not lead to severe violations of 
sum rules and can be used as a measure for the ghost 
problemJ^ Another, in our opinion, much more severe 
problem is the poor convergence at large energies. This 
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FIG. 3. (Color online) Comparison of the spectral weight 
calculated using the adaptive Lanczos-vector method with 
that calculated using the correction vector method for dif- 
ferent broadenings rj [U = 1.0, mcv = 300, tualm = 
3000, 250 Lanczos vectors) . 

could be addressed by spectral transformations, which 
are a standard method for improving convergence of ex- 
cited states within the Lanczos methodJ^ 

We now consider finite interaction, J7 = 1, and com- 
pare our method to the correction vector methodi^i^ The 
correction vector can be efficiently calculated by minimiz- 
ing a functional at every step of the DMRG,— but one 
must include a finite broadening in this process. Each 
broadening requires a separate DMRG run, in contrast 
to the ALM, which can evaluate the Green's function 
readily for any set of frequencies and arbitrary broaden- 
ing once the Lanczos coefficients have been calculated. 
In Fig. [3l we compare results from the two methods for 
two different values of the broadening, 77 = 0.1 and 0.2. 
We find very good agreement overall, especially for small 
energies. As expected, there are deviations at high en- 
ergies, which become more pronounced with decreasing 
broadening. 

We have shown that the ALM is capable of calculat- 
ing the spectral weights and poles accurately within the 



DMRG. In contrast to the OLM, we obtain the correct 
weights and poles of a Green's function up to energies 
of the order of half the bandwidth with only moderate 
resources, and, as long as the spectral weight is not too 
small, a reasonable reproduction of the spectrum even at 
larger energies. A clear advantage of the CFE in gen- 
eral is that it is based on an analytical expression for the 
Green's function, making it possible to evaluate it at an 
arbitrary set of frequencies with any possible broadening 
without having to repeat the DMRG calculation. In ad- 
dition, one can extract irreducible quantities such as self- 
energies directly as continued fractions. The inaccuracies 
that appear at higher energies can be traced to a loss of 
orthogonality of the Lanczos vectors due to systematic 
and numerical errors. In order to improve the accuracy of 
spectral functions at higher energies, the truncation num- 
ber m must be increased. This increase in m can become 
quite resource-intensive in standard DMRG implementa- 
tions. However, an improvement in efficiency could be 
achieved by using the Matrix-Product-State (MPS) for- 
mulation of the DMRG Within the MPS formulation, 
one can variationally optimize the three Lanczos vectors 
in the recursion formula separately, potentially leading 
to much better performance. Another advantage is that, 
within the MPS formulation, the Hamiltonian can be rep- 
resented exactly (see Ref.[l6|for an introduction), thereby 
further reducing systematic errors. Finally, we point out 
that this recursive approach to calculate a special basis 
for the evaluation of spectral functions is not limited to 
the Lanczos basis. For example, expansion in Chebyshev 
polynomials leads to a similar recursion formulaJ^ 

This work was supported by the DFG via SFB 602 and 
grant HO 2325/4-2. R.P. is supported by JSPS and the 
AvH. 

Note added: After completion of this work, we learned 
that Holzner et alr^ have shown that an adaptive method 
using Chebyshev polynomials in combination with MPS 
is highly efficient and gives accurate spectral functions. 
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